/*
 * DiffusionModel.java
 *
 * Copyright (C) 2002-2006 Alexei Drummond and Andrew Rambaut
 *
 * This file is part of BEAST.
 * See the NOTICE file distributed with this work for additional
 * information regarding copyright ownership and licensing.
 *
 * BEAST is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as
 * published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 *  BEAST is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with BEAST; if not, write to the
 * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
 * Boston, MA  02110-1301  USA
 */

package dr.evomodel.continuous;

import dr.inference.model.AbstractModel;
import dr.inference.model.Model;
import dr.inference.model.Parameter;
import dr.math.MathUtils;
import dr.xml.*;
import org.w3c.dom.Document;
import org.w3c.dom.Element;

/**
 * A class that can calculate the likelihood of a diffusion rate given a distance and time.
 * @author Alexei Drummond
 * @version $Id: DiffusionModel.java,v 1.5 2005/01/06 14:46:36 rambaut Exp $
 */

public class DiffusionModel extends AbstractModel {
	
	public static final String DIFFUSION_PROCESS = "diffusionProcess";
	public static final String DIFFUSION_CONSTANT = "D";
	public static final String BIAS = "mu";	
	
	/**
	 * Construct a diffusion model.
	 */
	public DiffusionModel(Parameter diffusionRateParameter) {
	
		super(DIFFUSION_PROCESS);
		
		this.diffusionRateParameter = diffusionRateParameter;
		addParameter(diffusionRateParameter);
		diffusionRateParameter.addBounds(new Parameter.DefaultBounds(Double.POSITIVE_INFINITY, 0.0, 1));
	}
	
	/**
	 * Construct a diffusion model.
	 */
	public DiffusionModel(Parameter diffusionRateParameter, Parameter biasParameter) {

		super(DIFFUSION_PROCESS);
		
		this.diffusionRateParameter = diffusionRateParameter;
		this.biasParameter = biasParameter;
		addParameter(diffusionRateParameter);
		addParameter(biasParameter);
		diffusionRateParameter.addBounds(new Parameter.DefaultBounds(Double.POSITIVE_INFINITY, 0.0, 1));
	}
	
	public double getD() {
		return diffusionRateParameter.getParameterValue(0);
	}
	
	/**
	 * @return the log likelihood of going from start to stop in the given time
	 */
	public double getLogLikelihood(double start, double stop, double time) {
		
		double D = diffusionRateParameter.getParameterValue(0);
		double bias = getBias();
		
		// expected variance of distances of given time
		double Dtime = D * time; 
		
		double unbiasedDistance = (stop - start) - (bias*time);
		
		
		//System.out.println("distance=" + unbiasedDistance + " time=" + time);
		
		// the log likelihood of travelling distance d, in time t given diffusion rate D
		return -0.5 * Math.log(Dtime) - ((unbiasedDistance*unbiasedDistance)/(Dtime));
	}
	
	/**
	 * simulate the diffusion process forward in time.
	 * @return the new value of the trait after given time.
	 */
	public double simulateForward(double value, double time) {
		double D = diffusionRateParameter.getParameterValue(0);
		double delta = MathUtils.nextGaussian();
		
		delta *= Math.sqrt(D * time);
		delta += getBias()*time;
	
		return value + delta;
	}
	
	
	/**
	 * @return the bias of this diffusion process.
	 */
	private double getBias() {
		if (biasParameter == null) return 0.0;
		return biasParameter.getParameterValue(0);
	}
	
	// *****************************************************************
	// Interface Model
	// *****************************************************************
		
	public void handleModelChangedEvent(Model model, Object object, int index) {
		// no intermediates need to be recalculated...
	}
	
	public void handleParameterChangedEvent(Parameter parameter, int index) {
		// no intermediates need to be recalculated...
	}
	
	protected void storeState() {} // no additional state needs storing
	protected void restoreState() {} // no additional state needs restoring	
	protected void acceptState() {} // no additional state needs accepting	
	protected void adoptState(Model source) {} // no additional state needs adopting	
	
	
    // **************************************************************
    // XMLElement IMPLEMENTATION
    // **************************************************************

	public Element createElement(Document document) {
		throw new RuntimeException("Not implemented!");	
	}
	
	// **************************************************************
    // XMLObjectParser
    // **************************************************************

	public static XMLObjectParser PARSER = new AbstractXMLObjectParser() {
		
		public String getParserName() { return DIFFUSION_PROCESS; }
			
		public Object parseXMLObject(XMLObject xo) throws XMLParseException {
			
			XMLObject cxo = (XMLObject)xo.getChild(DIFFUSION_CONSTANT);
			Parameter diffusionParam = (Parameter)cxo.getChild(Parameter.class);
				
			Parameter biasParam = null;
			if (xo.hasAttribute(BIAS)) {
				cxo = (XMLObject)xo.getChild(BIAS);
				biasParam = (Parameter)cxo.getChild(Parameter.class);
			}

			if (biasParam == null) {
				return new DiffusionModel(diffusionParam);
			}	
			return new DiffusionModel(diffusionParam, biasParam);
		}
		
		//************************************************************************
		// AbstractXMLObjectParser implementation
		//************************************************************************

		public String getParserDescription() {
			return "Describes a diffusion process.";
		}
		
		public XMLSyntaxRule[] getSyntaxRules() { return rules; }
		
		private XMLSyntaxRule[] rules = new XMLSyntaxRule[] {
			new ElementRule(DIFFUSION_CONSTANT, 
				new XMLSyntaxRule[] { new ElementRule(Parameter.class) }),
			new ElementRule(BIAS, 
				new XMLSyntaxRule[] { new ElementRule(Parameter.class) })
		};
	
		public Class getReturnType() { return DiffusionModel.class; }
	};
	
	// **************************************************************
    // Private instance variables
    // **************************************************************

	private Parameter diffusionRateParameter;
	private Parameter biasParameter;
}
